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ABSTRACT 

Under the assumption that LS 5039 is a system composed by a pulsar rotating around an 06.5V star in a 
~ 3.9 day orbit, we present the results of a theoretical modeling of the high energy phenomenology observed 
by the High Energy Stereoscopy Array (H.E.S.S.). This model (including detailed account of the system 
geometry, Klein-Nishina inverse Compton, 7-7 absorption, and cascading) is able to describe well the rich 
observed phenomenology found in the system at all timescales, both flux and spectrum- wise. 

Subject headings: X-ray binaries (individual LS 5039), 7-rays: observations, 7-rays: theory 



1. INTRODUCTION 

A few binaries have been identified as variable very-high- 
energy (VHE) 7-ray sources (Aharonian et al. 2005 a,b; 2006, 
Albert et. al. 2006, 2007). For LS 5039, a periodicity in the 
7-ray flux, consistent with the orbital timescale as determined 
by Casares et al. (2005), was found with amazing precision 
(Aharonian et al. 2006). Short timescale variability displayed 
on top of this periodic behavior, both in flux and spectrum, 
was also reported. Indeed, the H.E.S.S. observations of LS 
5039 (~ 70 hours distributed over many orbital cycles, Aha- 
ronian et al. 2006) constitute one of the most detailed datasets 
of high energy astrophysics. Accordingly, it has been sub- 
ject of intense theoretical studies (e.g., Bednarek 2006, 2007; 
Bosch-Ramon et al. 2005; Bottcher 2007; Bottcher & Demer 
2005; Dermer & Bottcher 2006; Dubus 2006a,b; Paredes et al. 
2006; Khangulyan et al. 2007). However, the rich observed 
phenomenology of the system has precluded an unifying ex- 
planation. 

The discovery of a jet-like radio structure in LS 5039 and 
the fact of it being the only radio/X-ray source co-localized 
with an EGRET detection, prompted a microquasar interpre- 
tation (Paredes et al. 2001). However, the current findings at 
radio and VHE 7-rays in the cases of LS I +61 303 (Dhawan 
et al. 2006, Albert et al. 2006) or PSR B 1259-63 (Aharonian 
et al. 2005), gave the perspective that all three systems are 
different realizations of the same scenario: a pulsar-massive 
star binary. Dubus (2006a,b) has studied these similarities. 
He provided simulations of the extended radio emission in the 
case of LS 5039, showing that the features found in high res- 
olution radio observations could also be interpreted as the re- 
sult of a pulsar wind. This result, the rest of similarities found 
with systems for which a pulsar companion is known, and the 
assessment of the low X-ray state (Martocchia et al. 2005), 
makes the pulsar hypothesis tenable. High energy emission 
from pulsar binaries have been subject of study for a long time 
(e.g., Moskalenko et al. 1993, Moskalenko 1995; Arons & 
Tavani 1993; Maraschi & Treves 1981; Romero et al. 2001). 
But the dichotomy between pulsars or black hole binaries is 
not settled (see, e.g., Mirabel 2006). Here, for testing whether 
the H.E.S.S. observations could be explained in a pulsar sce- 
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nario, we have assumed that the compact object in LS 5039 
is in fact a pulsar. Under this assumption, this Letter presents 
the results of a theoretical modeling of LS 5039 which is able 
to describe reasonably well the phenomenology found in the 
system. 

2. THE MODEL 

Following Sierpowska & Bednarek (2005) (see also 
Bednarek 1997), we consider that the volume of the binary 
system is separated by a shock wave into two regions with 
different properties, with the shock appearing as a result of 
the collision between the pulsar and the massive star winds. 
The position of the shock is defined by the parameter 77 = 
Lsb/(MV w c), where the wind velocity, V w , depends on the 
radial distance from the massive star, M is the star mass-loss 
rate, and Lsd is the spin-down luminosity of the pulsar (Ball 
& Dodd, 2001). Note that for 77 < 1 the star wind domi- 
nates over the pulsar's: the termination shock wraps around 
it. It is assumed that the pulsar and stellar winds are sym- 
metric. In the case of LS 5039 (see Table 1) the value of 77 
is between 0.5 (periastron) and 0.3 (apastron). Between the 
pulsar and the shock, inside the pulsar wind zone (PWZ), we 
assume that relativistic leptons are frozen in the magnetized 
pulsar wind which propagate radially from the pulsar. There- 
fore, their synchrotron losses are neglected in our cascade cal- 
culations. We consider then (Klein Nishina) inverse Compton 
(IC) cascades in this PWZ, assuming that the thermal radia- 
tion from the massive star dominates in this region (see e.g., 
Ball & Kirk 2000). Synchrotron emission are typically con- 
sidered in models where the 7-emission mainly comes from 
the shock region, especially applicable in the case of bina- 
ries with more eccentric orbits (e.g., Kirk, Ball & Skjaer- 
aasen 1999). In such very eccentric massive binaries with 
large separations the radiation field of the massive star does 
not dominate over the whole orbit and IC scattering is not ef- 
fective in the PWZ. We assume that the cascade initiated by 
a primary lepton in the PWZ develops in one dimension, i.e. 
in the direction of propagation of the primary particle. The 
charged products of this cascade arrive finally to the shock re- 
gion in the pulsar wind and follow the flow along the shock 
surface. Due to the magnitude of the opacity to IC and 7- 
7, most electrons effectively interact in the PWZ. Indeed, for 
instance for inclination i = 30°, for 1 TeV (100 GeV) elec- 
trons, the opacity to IC is ~ 1.5 (~ 9) in the PWZ. The en- 
ergy carried by electrons that actually reach the shock to be 
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trapped by the magnetic field is less than ~ 10% of the pri- 
mary injected energy. The secondary 7-rays move into the 
massive star wind region (MSWR). Some of them escape the 
binary system but a significant part can be absorbed due to 7- 
7 process. These cascades are followed by means of a Monte 
Carlo procedure (explained in general terms in Sierpowska & 
Bednarek 2005). The details of the LS 5039 implementation 
worked out for this paper (including the detailed account of 
the system geometry, the computation of opacities to the dif- 
ferent processes as a function of orbital phase, and many other 
details specific to our study) will be presented elsewhere. We 
assume that the pulsar-provided injected power in relativistic 
electrons is a fraction of the spin-down power, and that they 
are described by a power-law in energy. Table 1 presents the 
few free parameters of the model. 3 Additional model parame- 
ters are needed, but these are bound to fixed values (given also 
in Table 1) under multiwavelength observations (see, e.g., the 
papers quoted at the introduction). 

3. RESULTS 

Figure Q] shows the phase folded H.E.S.S. data for the inte- 
gral flux above 1 TeV (Aharonian et al. 2006) and the results 
of our theoretical model. Casares et al.'s (2005) ephemeris 
was used to fold the observational data. Each experimen- 
tal point represents a typical time span of observation of 28 
min. The corresponding phases for apastron, periastron, in- 
ferior (INFC), and superior (SUPC) conjunction are marked. 
Two periods and two different inclinations angles are shown. 
The trends in the H.E.S.S. datapoints are reproduced: orbital 
periodicity, a non-zero emission close to periastron, and the 
maximum of the fluxes is around INFC. 

VHE gamma-rays produced close to the star suffer se- 
vere absorption via pair production. The absorption process 
strongly depends on the direction of propagation with respect 
to the massive star. Indeed, the cross-section for pair produc- 
tion depends on the angle 6 between the VHE gamma-ray and 
optical photons, and has an energy threshold ex 1/(1 — cos 0). 
The level of absorption therefore depends on geometry and 
leads to orbital modulation of the VHE gamma-ray flux. For 
LS 5039, apastron and periastron are near INFC and SUPC. 
Absorption provides flux maxima at INFC, when 7-photons 
propagate towards the observer outward the system, and opac- 
ities are the lowest (cos — > 1). Flux minima happen at 
SUPC, when photons propagate toward the massive star (ab- 
sorption more effective, opacities the highest), as earlier dis- 
cussed, e.g., by Dubus (2006a). Absorption alone would how- 
ever produce strict modulation (zero flux) in the energy range 
0.2 to 2 TeV whereas the observations show that the flux 
at rsj 0.2 TeV is stable. Additional processes, i.e. cascad- 
ing, must be considered to explain the spectral modulation. 
Our model have these processes consistently included and the 
lightcurve details arise then as the interplay of the absorption 
of 7-rays with the cascading process, in the framework of a 
varying geometry along the orbit of the system. It can be seen 
that the differences between the two lightcurves correspond- 
ing to different inclinations are not large, but it is the detail at 

3 The first two free parameters shown are obviously related. However, 
note that the spin-down power defines the position of the shock in the system, 



INFC what is worth noticing. 

The spectra for LS 5039 was presented in two broad or- 
bital phase intervals around INFC and SUPC (Aharonian et 
al. 2006). A clear tendency for a harder spectrum at higher 
flux level was found. This is reproduced by the theoretical 
model herein presented. In order to compare with the broad 
phase-intervals spectra presented by H.E.S.S., we have also 
computed the phase-average of our predictions, see Figure 1. 
We found an overall agreement. H.E.S.S. observations have 
also provided the evolution of the normalization and slope of 
a power-law fit to the 0.2-5 TeV data in 0.1 phase-binning 
(Aharonian et al. 2006). The variability of the spectral index 
and normalization of these fits along the orbit is impressive. 
The use of a power-law fit has to do with low statistics in 
such shorter sub-orbital intervals: higher-order functional fit- 
tings (such as a power-law with exponential cutoff) provide a 
no better fit and were not justified. We can directly compare 
with these H.E.S.S. results fitting our spectral predictions in 
the same way. However, we note that particularly for some 
phases close to periastron, the theoretical spectrum is badly 
mimicked by a power-law (see the examples of Figure[2]). Us- 
ing these fits to each of the theoretically predicted spectra as a 
function of phase along the orbit, we compare with H.E.S.S. 
observations in Figure [3] The overall similarity of the the- 
oretical and observational fittings is notable. Depending on 
inclination there are a few points around INFC for which the 
observational data seem to be described by a harder spectrum 
than what we obtain as a result of fitting the theoretical pre- 
dictions. This can be modeled with a slightly harder injection 
spectrum for electrons. The normalization data point close to 
periastron is missed by our predicted fitting ranges. This is 
the result of our spectrum being badly fitted by a power-law 
at these phases. 



4. CONCLUDING REMARKS 

A detailed modeling of LS 5039, under a pulsar scenario, 
accounting for the system geometry, Klein-Nishina IC, 7- 
7 absorption, and a Monte Carlo computation of cascading, 
describes well the phenomenology found in the system at 
all timescales, both flux and spectrum- wise. The predictive 
power of the model is impressive: it is based on just a handful 
of free parameters. The results shown here were not obtained 
varying these parameters searching for a best fit, but rather 
exploring the predictions under minimal assumptions, what 
enhances their value. The overall agreement found does not 
however constitute a proof that the system is indeed formed 
with a pulsar companion (although contributes to circumstan- 
tial evidence in this direction). This agreement does provide 
a framework in which the pulsar assumption can be precisely 
tested with future observations. 

We acknowledge W. Bednarek, J. Cortina, I. Ribas, J. Rico, 
and N. Sidro for comments, extended use of IEEC-CSIC par- 
allel computers cluster, and support by grants AYA 2006- 
00530 and CSIC-PIE 2007501029. 

through 77, which also relates the stellar mass-loss rate M and Lsd- 
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TABLE 1 
Model Parameters 



Meaning 


Symbol 


Adoptee 


value 


Free model parameters 


Spin-down power of assumed pulsar 


^SD 


10 3Y erg s 


— i 


Fraction of spin-down power injected in relativistic electrons 


P 


10- 2 




Slope of the injected power-law description of the electron spectrum 


r 


-2.0 




VHE cutoff of the injection spectra (consistent with upper end of 7-spectrum in H.E.S.S. results) 




50TeV 




Inclination of the orbit of the system with respect to the line of sight (consistent with Casares et al. 2005) 


i 


30°, 60° 




Most important additional parameters, fixed by multiwavelength observations 


Radius of star 


R* 


9.3 R Q 




Mass of star 


M* 


23 M© 




Temperature of star 


T* 


3.9 x 10 4 


K 


Mass loss rate of star 


M 


10- 7 M 


yr _1 


Stellar wind (SW) termination velocity 




2400 km s 


-1 


SW initial velocity (SW velocity at radius r from the star: V w (r) = Vo + (Voo - Vb)(l - R+fr) 1 - 5 ) 


Vo 


4kms _1 




Distance to the system 


D 


2.5 kpc 




Eccentricity of the orbit 


8 


0.35 




Semimajor axis 


a 


0.15 AU - 


- 3.5 ft* 


Longitude of periastron 


Up 


226° 






Fig. 1. — Left: H.E.S.S. run-by-run folded observations of LS 5039 with the results of the theoretical model for two different inclination angles. Right: 
H.E.S.S. VHE spectra around INFC and SUPC together with theoretical predictions in equal phase intervals. 



Albert J. et al. 2007, ApJ Letters 665, 51 
Arons J., & Tavani M. 1993, ApJ, 403, 249 
Bednarek W. 1997, A&A 322, 523 
Bednarek W. 2006, MNRAS 368, 579 
Bednarek W. 2007, astro-ph/06 11291 
Ball, L. & Dodd, J. 2001, PASA, 18, 98 

Bosch-Ramon, V., Paredes, J. M., Ribo, M., et al. 2005, ApJ, 628, 388 
Bottcher, M. 2007, Astroparticle Physics 27, 278 
Bottcher, M., & Dermer, C. D. 2005, ApJ, 634, L81 

Dhawan V., Mioduszewski A. & Rupen M. 2006, Proceedings of Science 
(MQW06) 052, VI Microquasar Workshop, Sept. 18-22, 2006, Como, 
Italy. 

Dermer, C. D. & Bottcher, M. 2006, ApJ 644, 409 

Casares, J., Ribo, M., Ribas, L, et al. 2005, MNRAS, 364, 899 

Dubus G. 2006a, A&A 451, 9 



Dubus G. 2006b, A&A 456, 801 

Khangulyan D., Aharonian F. A., Bosch-Ramon V. 2007, arXiv:0707.1689 

Kirk J. G., & Ball L., Astropart. Phys. 2000, 12, 335 

Kirk J. G., Ball L. & Skjaeraasen O., Astropart. Phys 1999, 10, 31 

Maraschi L., & Treves A. 1981, MNRAS, 194, 1 

Martocchia, A., Motch, C., Negueruela, I. 2005, A&A 430, 245 

Moskalenko, I.V, Karakula, S., Tkaczyk, W. 1993, MNRAS 260, 681 

Moskalenko, I.V. 1995, Space Science Rev. 72, 593 

Paredes, J. M., Marti, J., Ribo, M., & Massi, M. 2000, Science, 288, 2340 

Paredes, J. M., Bosch-Ramon, V, Romero, G. E. 2006, A&A 451, 259 

Romero G. E, Kaufman Bernado M. M., Combi J. A. & Torres D. F. 2001, 

A&A 376, 599 
Sierpowska A. & Bednarek W. 2005, MNRAS 356, 711 



4 



A. Sierpowska-Bartosik & D. F. Torres 




Fig. 2. — Examples of the theoretical results, for different inclinations, for the spectrum of single phases fitted with a power law in the range 0.2-5 TeV. 





Fig. 3. — Parameters of power-law fittings to the spectra as a function of phase both for the observational H.E.S.S data and our theoretical curves. We show 
two sets of results corresponding to fitting in different energy ranges (0.2-1 TeV, 0.2-5 TeV), and inclination angles (top: i = 60°, bottom: i = 30°). The size 
of the shading represent the error in the parameter fitting of our model. 



